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Abstract 



We study four dimensional SU(2) lattice gauge theory in the Hamiltonian 
formalism by Green's Function Monte Carlo methods. A trial ground state 
wave function is introduced to improve the configuration sampling and we 
discuss the interplay between its complexity and the simulation systematic 
errors. As a case study, we compare the leading strong coupling approximation 
and an improved 4 parameters wave function with lxl and 1x2 plaquette 
terms. Our numerical results favors the second option. 

PACS numbers: 11.15.Ha, 12.38.Gc 
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I. INTRODUCTION 



In the study of Lattice Gauge Theory, a major alternative to the usual Lagrangian 
formulation is the Hamiltonian version of Kogut-Susskind 0. This approach is relevant 
both for theoretical reasons, e.g. universality checks, and for practical purposes because it 
seems the natural choice for static problems, like spectroscopy 0. 

An important feature of the Hamiltonian formulation is that it can exploit analytical 
approximations of the ground state wave function to improve the numerical evaluation of 
physical quantities 0. This well known procedure is called Importance Sampling J4|. How- 
ever, even if an accurate wave function is certainly welcome, on the other hand it can also 
require expensive calculations. A careful analysis is therefore necessary to understand the 
interplay between accuracy and performance and to determine which is the optimal com- 
plexity of the trial wave function. 

A recent investigation of Hamiltonian methods in the numerical study of non Abelian 
gauge theories can be found in |J where the SU(2) model is studied in 2+1 dimensions. Here, 
we focus mainly on the wave function issue and study the 3 + 1 dimensional SU(2) pure gauge 
lattice theory by a particular Green's Function Monte Carlo (GFMC) || algorithm that 
computes statistical averages over an ensemble with a fixed number of random walkers 
This algorithm has the desirable feature of allowing a particularly simple analysis of the 
systematic errors. 

On general grounds and motivated by strong coupling expansions j§| we consider trial 
wave functions which are expressed by combinations of gauge invariant Wilson loops with 
unknown coefficients a. Two main issues are addressed: the cost of computing large loops 
and the practical simultaneous optimization of many parameters a. About the latter, we test 
a recently proposed algorithm for the adaptive optimization of a. It is a method that has 
been successfully applied to the study of a very simple system with U(l) gauge symmetry || 
and that deserves a more detailed analysis on a realistic test-bed like the non Abelian SU(2) 
four dimensional model. 

As we shall see, the most serious systematic error is related to the finite population size, 
namely the number of walkers. The extrapolation to an infinite population is difficult and 
much easier when an improved wave function is used. Working with a four parameter trial 
function, we present results for the extrapolated quantities and study the scaling of the 
string tension comparing our results with similar Lagrangian calculations. 

The plan of the paper is the following: in Sects. f|, HU an< ^ 0' we P resen t the GFMC al- 
gorithm in its general form for SU(N) lattice gauge theories with Importance Sampling, fixed 
number of random walkers and adaptive optimization of the trial wave function. In Sec. |V], 
we apply the proposed algorithm to the SU(2) model in four dimensions and discuss the 
numerical results. 

II. REVIEW OF GFMC FOR SU(N) GAUGE THEORY 

In this Section, we recall the Feynman-Kac- Nelson formula for the matrix elements of the 
evolution operator associated with the SU(N) Kogut-Susskind Hamiltonian. This formula 
provides probabilistic representations of several interesting quantities. The GFMC algorithm 
is just one of its numerical implementations. 
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Let U = {U{\ be the set of group elements, one for each link /; let E = {Ef} be the 
associated electric field defined by 

[E?,U v ]=T a U l 5 ltV , (2.1) 

where {T a }, a = 1, . . . , N 2 — 1 are SU(N) generators. Finally, let V(U) be a gauge invariant 
real potential. The Kogut-Susskind Hamiltonian for SU(N) lattice gauge theory is then 
Hks — 9 2 H where g is the gauge coupling and 

1 1 N2 ~ l 

H = H + V(U), H = -E 2 = -Y1 E(^) 2 - (2-2) 

A / a=l I 

For the Wilson action, V(U) is proportional to the sum of the smallest lxl plaquette loops 

y(U) = ~^J2^(U p + Ul), (2.3) 

but the present discussion holds for a general V{U). 

The matrix elements of the Euclidean evolution operator exp(— tH) in the basis of Ui 
eigenstates admit the small t expansion 

D(U", U', e) = (U"\e- £H \U'} = D (U" } U', e)e- £V{ul) + 0(e 2 ) } (2.4) 

where the factor 

D (U",U',e) = (U"\e- £Ho \U f ), (2.5) 

can be interpreted as a probability density for the random transition U' — > U" due to the 
relations: 

D (U", U', e)>0 and J D (U", U', e)dU" = 1, (2.6) 

(dU is the SU(N) Haar measure). Inserting intermediate states between U" and U', we 
obtain as usual 

D{U", U',t) = J (U n \e-W N+1 » H \U N ) ■ ■ ■ (U^e-WW+WWdUi ■ --dU N , (2.7) 

and taking the N —* oo limit we recover the Feynman-Kac-Nelson path integral representa- 
tion oiD(U",U',t): 

D(U",U',t)= I VU(t)e-fo dTV W\ (2.8) 

JU(0)=U', U(t)=U" 

where the measure T>U(t) stands for the N — >• oo limit of the average Eq. ( |2.7| ) over random 
group polygonals (Uq — U',Ux, . . . , £/jy, Z7jv+i = U"). 

The potential V(U) fluctuates along the path U(t) and determines the error on the 
estimates of physical quantities. To see how this happens, let us consider in this framework 
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the calculation of the ground state energy E . The simplest procedure is to compute the 
limit 

E = lim E(t,U'), (2.9) 

where 

E(t, U') = -j t log J D(U", U\ t)dU", U' arbitrary. (2.10) 

It is easy to check that 

E(t, U>) = J JM=E U 1 UJ ee (VW))), (2-11) 

f u{0)=u/ VU(t)e-fo d ^)) 

where (•) is the average over free trajectories on the group manifold weighted by the expo- 
nential term. Hence, the fluctuations of V control the noise in any estimator of Eq. Other 
problems in the actual numerical evaluation of (V) can be discussed separately and will be 
considered in the next Section. 

The introduction of a trial wave function is a clever way to reduce the fluctuations of 
V. To this aim, one considers a positive function ^o(U) = expF(U) and the unitarily 
equivalent Hamiltonian 

H = *o[H + V(U)]*o 1 , (2.12) 

which turns out to be 

H = Ho + V, (2.13a) 
Ho = -E 2 + iE-VF, (2.13b) 
E-VF^Y.EfV^F (2.13c) 

l,a 

y = y-\ E(v«^) 2 -\e v^f, (2.13d) 

a, l a, I 

where V Qj j is the group invariant derivative on SU(N) associated with link /. The following 
approximate factorization of the H propagator holds 

D = D e- £V + 0{e 2 ), (2.14) 
where D is the kernel of exp(— eH ) and satisfies the fundamental relations 

D (U",U',e) > and J D (U", U', e)dU" = 1, (2.15) 

(the important point is that H is normal ordered with all the E operators acting from the 
left). As in the F = case, the finite time propagator D may be expressed in terms of a 
weighted average 
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D(U",U't)=[ VUit) e -fo drV W\ (2.16) 

JU(p)=U', U(t)=U" 

where VU is the measure determined by the form of D at infinitesimal times. The explicit 
first-order discrete construction of the paths is done by choosing a time step e and updating 
U n — ► U n+ i according to the rule 

U n+1 = U R {e) ■ U D (s, U n ) ■ U n , (2.17) 

where Ur(e) is a random SU(N) element distributed according to the heat kernel |]rO| , |TT] at 
time e and Ud{e, U) is the drift 

U d (e, U) = exp(ieT ■ VF(U)). (2.18) 

When happens to be an exact eigenstate of H with eigenvalue E, we have V = E 
and the formula 

E(t) = (V(U(t))) = E, (2.19) 

provides the correct eigenvalue with zero variance, namely no statistical error. In other 
words, if we require V to be a constant, then the problem is completely equivalent to solve 
the Schrodinger equation for the Kogut-Susskind Hamiltonian. In the following, we shall 
try to obtain a constant V at least in some approximate sense. 

In the above discussion we focused on the ground state energy which is the simplest 
observable. A more general case is that of ground state matrix elements of diagonal operators 
in the basis of the link variables U : 

(Q) = (0\n(U)\0), (2.20) 

where |0) denotes the ground state of H. They can be evaluated by computing 

(O) = lim lim £l(t 2 ,t 1 ,U'), (2.21) 

t<2 — > + 00 t\ — > + 0O 

where 

r*i+*2 



S u(fi)=u ,VU{t)e-K a *^MV(tf(tO) 



n(*2, h, u') = ^ ftl+ « a , • ( 2 - 22 ) 



I m=u ,VU(t)e-Io 



drV{U(r)) 



A similar formula can be also devised for many time correlation functions of operators in 
the Heisenberg representation. 



III. STOCHASTIC RECONFIGURATION 

In this Section, we discuss the actual calculation of the above mentioned averages over 
trajectories. Here, the trial wave function wm n °t play any role. 

A naive implementation of the Feynman-Kac- Nelson representation generates group tra- 
jectories according to the measure T>U and weights them with a potential dependent weight. 
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Such weight increases or vanishes exponentially fast with respect to the trajectory length t. 
This is a serious numerical problem referred to as the "exploding variance problem". The 
name stems from the fact that the weight variance of a collection of random walkers explodes 



as t — > +00 overwhelming the numerical capabilities of any machine jn| . 

A possible solution to this problem consists in killing and cloning the random walkers 
with definite rules in order to delete the walkers with low weights and duplicate the others. 
The main disadvantage is that the size of the walker population varies in time and may 
become unstable. 

The recently proposed Stochastic Reconfiguration Algorithm 0] is a simple implemen- 
tation of the killing and branching idea, but with the desirable feature of dealing with a 
population with a fixed number of walkers. Its coding is therefore much simpler. 

To this aim, a finite collection of K walkers, an ensemble, is introduced 

£ = {(U^ n \t),^ n \t))} 1 < n < K , (3.1) 
where the weights are defined by 

u {n) {t) =exp- f dTV(U {n \r)). (3.2) 
Jo 

The variance of the weights over the ensemble £ is 

W(t)=Varuj(t) = ±j:(J k \t)f-(±f:J k \t)) . (3.3) 

k—1 \ k=l I 

The average of a function f(U) over £ is defined as 

(f)e = ^j:f(U ik) W k \ (3.4) 

k=l 

When W(t) becomes too large, our aim is to transform £ in a new ensemble £' with a rule 
£—>■£' such that the variance W is reduced while the averages are kept constant. Actually, 
this is a formidable task, but we can achieve it at least approximately, namely in the K — > 00 
limit. To this aim, we extract K independent integers n>i, . . . , n K with n» G {1, . . . , K} and 
probabilities 

Prob(n; = k) = — ^ -g-. (3.5) 



We then define £' by 



£' = {(U^\t),l)} 1 < l < K . (3.6) 



In other words, we replace the information carried by the weights with that encoded in the 
multiplicities of unit weight walkers. By the large numbers law we expect 

(f)s - (f)e> A -=5° (3.7) 
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In actual calculations, there is therefore a systematic error which is removed by performing 
simulations with different values of K and somehow extrapolating the K — > oo limit. 

The calculation of observables when stochastic reconfiguration is in progress does not 
present any additional problem. The only point we stress is that, for a general operator Q, 
the weights at time t\ + t 2 are used to weight Vt{t{) and therefore one must keep track for 
each reconfigured walker of its states at that past time (usually, a small number of updates 
is enough). 



IV. ADAPTIVE OPTIMIZATION OF THE TRIAL WAVE FUNCTION 

Let us consider a trial wave function ^/q(U,sl) depending on some parameters 
) and the GFMC calculation of an observable quantity that, for simplicity, 
we choose to be the ground state energy E . 

We do not discuss the systematic error associated with e appearing in Eq. ( |2.17| ). The 
choice of e is related to the approximation involved in the operator splitting and also in 
the evaluation of the SU(N) heat kernel. A common model independent choice is to take e 
enough small to approximate the diffusion over SU(N) with that on the flat tangent space at 
the identity. Problems may then occur only if V-F becomes too large, a constraint that can 
be checked during the simulation. In the following we shall choose a reasonably small e and 
work under the hypothesis that the e — > extrapolation does not change our conclusions on 
the trial wave function dependence. 

A simulation with a population of K walkers will provide after S Monte Carlo steps an 
approximate estimator E (S, K, a) of £"0- This is a random variable with the property 

(Eo(S,K,a)) = Eo + CLiK,*), (4.1) 

and 

VarE (5,K,a) = ^^, (4.2) 

where (•) denote the average over the Monte Carlo realizations. The finite population 
functions Ci(K, a) satisfy 

lim Ci(K,a) = 0, (4.3) 

where, by a self- averaging argument, it is reasonable to guess the asymptotic form C2{K, a) ~ 
j(a)/y/K for large K. The average of E extrapolated at K — ► 00 is exact and indepen- 
dent on the trial parameters a. On the other hand, the function C2(K, a) may be strongly 
dependent on them. Since c 2 determines the statistical error at fixed computational time, it 
is important to minimize it. 

If the family of functions ^o(U, a) includes the exact ground state at the special point 
a = a* then we know that 

d(K,a.*) = c 2 (fr,a*) =0. (4.4) 

In a less optimal situation one has the two possibilities of minimizing c\ or c%. Motivated 
by [13|], we choose to pursue the second alternative and seek a minimum of c<i- 
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Moreover, we propose to optimize a adaptively. This can be done in the following way: 
the constant C2 is related to the non zero fluctuations of V defined by 

V(U, a) = V(U) - \{VF{U, a)) 2 - l -V 2 F{U, a), (4.5) 

where ^o(U, a) = exp F(U, a). The ensemble variance of V is a function of a 

k / 1 K 



Var £ V = ± £ V(U {h \ *?^ k) -(jrT, V(U {k \ *)u<A . (4.6) 

k=l \ k=l J 



Our proposal is to update a according to the simple equation 

a n+ i = a„ - r]V a Var£- V, (4.7) 
where rj is a constant parameter and a n is the value of a at the n-th update. In other words, 



following the spirit of ||13| , we implement a local minimization of the weight variance as a 
driving mechanism for the free parameters. It is important to remark that it is incorrect 
to look for the special parameters a that minimize Var^ V for a fixed ensemble £. In fact, 
the equilibrium distribution of the link variables in the ensemble do depend on a and the 
addition of Eq. ( |4.7| ) to the Monte Carlo core establishes a non trivial feedback between the 
two sets of variables. 

The coupled set of equations Eqs. Q2.17| , f4.7| ) for the evolution of a and the random walkers 
is non linear and discrete. The above procedure is reasonable, but a test in explicit examples 
is required to check the convergence and stability of the method. A first analysis is described 
in for a toy model with U(l) gauge symmetry in low dimensions. Here we consider the 
less trivial non Abelian four dimensional case. 



V. SIMULATION AND RESULTS 
A. General Setup 

We simulate the SU(2) lattice gauge theory on a spatial 4 3 lattice. The temporal ex- 
tension is in principle infinite. Of course, the meaning of this statement is that in the 
Hamiltonian formulation we are projecting onto the exact ground state since we are ap- 
plying many times an approximate form of exp(— eH). In all our simulations e = 0.01. 
Another parameter that we keep fixed for simplicity is the number of elementary time evo- 
lutions between two stochastic reconfigurations. We choose 10 steps which is roughly the 
decorrelation time of the lattice configuration. A priori, a better choice of this parameter 
can improve the convergence of the K — > 00 extrapolation, but we checked that in our range 
of parameters, the sensitivity on the choice of K is quite small. We also fix the learning 
parameter 77 = 0.001 for the adaptive determination of the wave function coefficients. 

The parameters that remain free are therefore: the number K of walkers, the number of 
coefficients in the ground state wave function and of course the coupling constant g which is 
needed to discuss the possible scaling of the measured quantities. In more details, for each 
value K = 2, 10, 50, 100 and for each choice of the trial wave function we run the algorithm 
at five different coupling constants 



S 



A = A/g 4 = 0.8, 0.9, 1.0, 1.1, 1.2, 



(5.1) 



and measure the ground state energy E and all the Wilson loops Wjj with area I J < 9. 
We recall that the variable 4/g 4 is natural in the strong coupling expansion analysis. 



B. Adaptive Optimization 

The first topic that we discuss is the adaptive optimization of the free parameters a in 
the trial wave function. Motivated by the long wavelength strong coupling analysis of ||, 
we consider the following form of F 

F = c x W Xxl + c 2 W (lxl) 2 + c 3 W 1 

X 2, planar X 2, bent; 

(5.2) 

where Wi x i, W^ixi) 2 , W / ix2, P ianar and VFi^.bent are the sum of all the plaquettes of the kind 
shown in Fig. [TJ The coefficients ci, C2, C3 and C4 may be approximated as in || with a 
power series in 1/g 2 computed on the infinite lattice. Here, we exploit the dynamical sys- 
tem Eq. ([4.7|) to determine the finite size optimal Cj. In the analysis we discuss three specific 
scenarios that we characterize by the level / = 0, 1 or 4, the number of free parameters: 

1 = 0: F = 0, no Importance Sampling; 

1 = 1: c 2 = c 3 = c 4 = 0, one free parameter; 

I = 4: all Q are optimized. 

In Fig. we show for example the evolution of C\ with Monte Carlo time in the case 
K = 2, 1/g 2 = 0.5 and 1 = 1. The coefficients q typically stabilize on values adaptively 
computed by the algorithm with small fluctuations that can be averaged in the long time 
limit. The initial values of c« for K = 10, 50, 100 are taken from the run with the nearest 
smaller K in order to reduce thermalization times. The procedure turns out to be stable at 
each considered K and g. The fluctuations are bigger at smaller g, that is in the critical limit, 
while the equilibrium values agree reasonably well with the strong coupling approximation 
at large g. If the fluctuations were too noisy, then Stochastic Gradient Approximation 
techniques fTM could be applied by gradually reducing r\ in Eq. ([4.7|). 



C. K — > 00 extrapolation 

For each level /, we consider simulations with K = 2, 10, 50 and 100. The K — > 00 limit 
does not depend on /, but the rapidity of the convergence yes. In Figs. (|3|-[4|) we plot the K 
dependent estimates of E and Wi x i for / = 0, 1, 4 at 1/g 2 = 0.5. In the left plots we see 
clearly that for / = and I = 1 the extrapolation K — > 00 is largely uncertain. Unless we 
know the asymptotic behavior in K, we cannot safely determine the correct limit with the 
present data. A non rigorous guess is a power law behavior with finite K corrections of the 
form 

c i(^ a ) = fS + °(^" a )> a(a)>0, (5.3) 
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that has been observed empirically ||. Such behavior can be matched to the K = 10, 50, 100 
data, but its choice is rather arbitrary and, to be supported, would require more points at 
larger K. On the other hand, data computed in the I = 4 scheme show a plateau beginning 
at K ~ 50 and data at K = 100 will therefore be used as an estimate of the K — > oo limit. 

To check consistency of the runs at different K and to evaluate quantitatively the im- 
provement gained with larger K, we introduce an effective K 

K' = K/X(l), (5.4) 



depending on the rescaling factor A(/). In the right plots of Figs. (043); we show how the 
curves for Eq and Wi x i measured at different K collapse on a single curve when shown as 
functions of K' . The values of A(Z) are 

A(0) = 1.6- 10 3 , A(2) = 5, A(4) = l. (5.5) 

The raw simulation, I = 0, is very poor, as expected. On the other hand, the one parameter 
wave function gives the same accuracy of the / = 4 one when the number of random walkers 
is roughly 5 times bigger. 

We now determine the computational cost as a function of I. Let n be the number of 
Monte Carlo iterations; the statistical error in the measure of a given observable and the 
CPU time are respectively given by 

+ - nK a\ 

£stat — T= i t CPU — , {O-O) 

VnK Vi 

where 7; and v\ are normalization constants expected to be a decreasing function of /. Of 
course, 7; depends on the chosen observable. From Eq. (|5.6|) we deduce that 



7, 2 1 

topu = -— (5.7) 

V l £ stat 

is independent on K. From our data, we can give the rough estimate of the ratios V\/v± and 
71/74 when the coupling is 1/g 2 ~ 0.5. We find 

Vi/v A ~ 8, 71/74 ~ 2, (5.8) 

and therefore t C pu,i=i — 2 t C pu,i=i- Of course one must also remember that, as discussed 
above, I = 1 does not show a clear plateau when K < 100 and that very large K cannot be 
used in 3 + 1 dimensions because of storage limitations. Larger Wilson loops show a similar 
behavior with slightly better improvement in the case I = 4. 

In this particular model and at the considered volume and couplings we summarize the 
comparison between / = 4 and I = 1 at fixed systematic and statistical errors by saying 
that the improved wave functions are twice slower, but require a storage ~ 5 times smaller. 
In the following, we shall pursue the I = 4 choice because only in that case, our data can be 
safely considered asymptotic with respect to the K — ► 00 limit. 
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D. Scaling of the string tension 



The I — 4, K — 100 measurements of E , the Wilson loops with area < 9 and the 
coefficients c« are shown in Figs. (0-0) as functions of 1/g 2 . The energy measurements agree 
well with the second order strong coupling expansion 

E ° -7 ~ 5? + ° (£) (8 ' 9) 

An important check of the simulation requires an attempt to investigate scaling. The easiest 
observable quantity with simple Renormalization Group behavior is the string tension a that 
we obtain from Creutz's ratios 

2 W ItJ W I+ltJ+1 

a = hm a ItJ , a a I:J = - log— — . (5.10) 

i>J-*°° W I)J+ iW I+ ij 

This must be compared with the asymptotic scaling prediction 

1 /24vr 2 \ 51/121 / 12tt 2 \ 

Even if we are work on a lattice with rather small spatial extension, it is interesting to 
compare our results with those obtained in the Lagrangian formulation. We choose the 
measurements in [pl| as representatives of what can be obtained at values of the coupling 
in the range Eq. ( |5.1|) . In Fig. § we show a^fa\^ and f(g). Loops larger than 2x3 cannot 
be used since in our computation they have too large errors. 

The Lagrangian measurement is At = (1.19 ± 0.15) • 10~ 2 y/a. Taking into account the 
SU(2) A parameter ratio A^/A^ = 0.84 |16[, we obtain the band shown in the Figure. As 
one can see, the Hamiltonian measurement is slightly larger, but reasonably consistent with 
the Lagrangian one. 



VI. CONCLUSIONS 

The aim of this paper has been a study of the non perturbative behavior of the SU(2) 
lattice gauge theory in 3 + 1 dimensions by Hamiltonian Monte Carlo methods. Our results 
can be summarized into two main statements. 

First, we have shown that the algorithm || is able to adaptively optimize the many 
parameters ground state trial wave functions needed to guide GFMC simulations of realistic 
models. 

Second, we have discussed to what extent the use of improved wave functions is actually 
necessary, at least with the considered volume and couplings. From the purely computa- 
tional point of view our data show that wave function improvement reduces the algorithm 
performance, but permits a substantial gain in storage. This means that we are allowed to 
work with a smaller number of walkers, K. This can be an important advantage, especially 
on larger lattices or at smaller couplings where the minimum K at which a safe extrapolation 
to the K — > oo limit is feasible is expected to increase. 
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As a non trivial check of the algorithm and of the strategy for the removal of the finite 
K systematic error, we studied the string tension with a four parameter trial wave function 
and found a scaling behavior compatible with similar Lagrangian measurements. 

Our investigation has been performed on a cluster of small computers and we could not 
attempt to measure glueball mass ratios. At this stage, the Hamiltonian formulation does 
not offer any special advantage with respect to the Lagrangian one. The storage gain as- 
sociated to the absence of the temporal dimension is compensated by the memory required 
to manage the walkers ensemble. In practice, we did not examine finite size scaling of the 
algorithm, a task that we leave for future studies. We stress however that many of the 
advanced techniques that are currently used in Lagrangian calculations [0, can be applied 
in the Hamiltonian formalism as well. Examples are the use of refined glueball operators or 
the improvement programme that has been recently extended to Kogut-Susskind Hamilto- 



nians [|17 |. 
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FIGURES 



FIG. 1. Plaquettes appearing in the I = 4 trial wave function. A sum over the lattice is 
understood. 

FIG. 2. Monte Carlo history of the trial wave function coefficient c\ at 1/g 2 = 0.5 and K = 2. 
The initial state of the system is cold. 

FIG. 3. Dependence of Eq on K with and without reparametrization K — > K' for the three 
schemes I = 0, I = 1 and 1 = 4. The coupling constant is 1/g 2 = 0.5. 

FIG. 4. As Fig. g, but for the 1 x 1 Wilson loop. 

FIG. 5. Coupling constant dependence of Eq at / = 4, K = 100. 

FIG. 6. Coupling constant dependence of the trial wave function coefficients at I = 4, K = 100. 

FIG. 7. Coupling constant dependence of the Wilson loops at / = 4, K = 100. 

FIG. 8. Scaling plot of the string tension. Comparison between a^Ja\^ and the asymptotic 
Renormalization Group prediction. The band corresponds to the Lagrangian measurement as 
explained in the text 
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